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Abstract 

A detailed simulation of vertical showers in atmosphere produced by primary gam- 
mas and protons, in the energy range 1-100 TeV, has been performed by means of the 
FLUKA Monte Carlo code, with the aim of studying the time structure of the shower 
front at different detector heights. It turns out that the time delay distribution can be 
fitted using few parameters coincident with the distribution central moments. Such 
parameters exhibit a smooth behaviour as a function of energy. These results can 
be used both for detector design and for the interpretation of the existing measure- 
ments. Differences in the time structure between gamma and proton induced showers 
are found and explained in terms of the non-relativistic component of extensive air 
showers. 
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1 Introduction 



Recent papers concerning calculations of Extensive Air Showers (EAS)|M include results 
on time delay of particles from the shower front, in view of the comparison with existing 
recent measurements 0, [3L 0. As a matter of fact, the first measurements of the structure 
of the EAS front were attempted by Bassi, Clark and Rossi in 1953Q, but many other 
have contributed As far as the many simulation works on EAS are concerned, a 
part from the quoted ref. M , only a limited fraction of them has considered the question 
of the time structure filf]- [p3[| . 

The interest in this topic has been renewed by recent experimental data concerning 
the detection of an anomalous delayed component^]. In this framework it is therefore 
important to achieve a better and more detailed knowledge of this time structure. Fur- 
thermore, eventual fast simulation tools based of simple parametrisations constructed on 
this knowledge can be helpful not only for analysis, but also for design of new detectors, 
trigger evaluation, et cetera. 

When a full shower simulation is considered, the structure of the time front depends to 
a large extent on the details of particle transport algorithms. Therefore, it can be useful 
to consider the outputs of very detailed codes, and compare it with the results of the codes 
optimized for the cosmic ray physics, which often contain simplifications in order to reduce 
the computer time. A valuable study of the time structure of Extensive Air Showers, for 
very high energies, based on the use of the CORSIKA code[25] has been presented in [||. 

Here we present the study of the time structure of the e.m. secondary component 
generated by photon or proton, as calculated mainly with the FLUKA96 code. We have 
already used such code in pq| , where we published the parametrisations of e.m. size 
resulting from photon and electron induced sub-showers. Such code has already been 
successfully used in other cosmic ray applications |2^, Bq| . 

We have explored primary energy from 100 GeV up to 100 TeV, and injected at the 
top of the atmosphere. This energy range is of particular interest for gamma-astronomy 
purposes, mainly in the chance of possible differences in the time structure between gamma 
and proton initiated showers. 

After describing the MC set-up, we discuss the search for the best functional shape to 
fit the time distribution, and then results are discussed. Numerical results of characteristic 
parameters are given for different energies, secondary particles and detector altitudes. 



2 The MC codes 

Among the shower programs available in the High Energy Physics community, we have 
mainly considered FLUKA96 [p9| , although in some cases we have also used GEANT ver- 
sion 3.21 [pQ l, using the interface to the FLUKA package (1992 version) for hadronic 
interactions. It must be stressed that GEANT-FLUKA and FLUKA96 are not the same 
thing, for different reasons: the e.m. simulation code is different pl|], particle transport 
algorithms are different and, most of all, the FLUKA interface of GEANT contains only 
part (updated to 1993) of the hadronic interaction model of full FLUKA96. 

In particular, the full FLUKA96 code employs transport algorithms with refined path 
length corrections^] associated to multiple scattering, which are essential in problems 
involving low energy electrons and positrons. Furthermore, FLUKA96 is completely in 
double precision, thus allowing a much more accurate definition of a finely segmented 
geometry set-up which extends itself up to many tens of kilometres. Another advantage 
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coming from double precision is the timing accuracy even below 1 ns. An additional 
advantage of the full FLUKA96 package is the inclusion of very detailed models of nuclear 
excitation and low energy neutron transport. As we shall see in the following, this last 
topic has some importance for the argument of this paper. The two different codes adopted 
in this work provide slightly different results, in particular for logi^t /ns < 0, because of 
the insufficient precision of the GEANT transport algorithms. However, for logiot/ns 
greater than that, the differences in the bulk of arrival time distributions are substantially 
within statistical fluctuations. For all these reasons, in the following sections we shall limit 
ourselves to the discussion of the results obtained with FLUKA96. 

The atmosphere is defined by a stack of box volumes of rectangular basis and thickness 
increasing with the height above the sea level. Any volume corresponds to a depth of ~ 
24 g/cm 2 . In each box the density is uniform, and it is chosen in such a way that an 
approximation to the standard U.S. atmosphere is performed according to the Shibata 
fit|33|]. The chosen depth granularity in our approximation is about one half of radiation 
length (37.66 g/cm 2 ) in air. We have limited the top of atmosphere at 1 g/cm 2 , and the 
bottom is at 1025 g/cm 2 . Particles have been injected at the altitude H = 45.445 km 
above the sea level. We have used the same elemental composition at all depths. The 
kinetic energy cut for secondary charged particles has been fixed to to 1 MeV, while for 
secondary photons we have chosen 0.5 MeV in order to include the contribution from e + e~ 
annihilation. Particles have been recorded at three different detector heights, correspond- 
ing respectively at sea level, 1000 and 2000 m above sea level. These altitude values have 
been chosen in order to apply our work mainly to the case of the Gran Sasso laboratory 
site. For our investigation, we generated only vertical showers, at different log-spaced 
energies: 100, 177, 316, 562, 1000, 1778, 3128, 5623, 10000, 17780, 31620, 56230, 100000 
GeV. 



3 Analysis 

We present results for different groups of secondary particles. We have considered sepa- 
rately photons, electrons ( both e + and e~ ) and other charged particles (muon, pions, 
etc.). The last group of secondary is relevant only in the case of primary hadrons. It 
must be noticed in fact that the adopted codes did not activate the hadro-production by 
photons. 

As already done by other authors [Q] we defined as relevant variable, the time delay 
t with respect to a spherical front moving with light speed c, originating from a fixed 
injection point the atmosphere, i.e. we are giving the shower disk deviation from such a 
spherical front. The exact definition of our delay variable t is given by: 

t — t arrival ~ (1) 

where r is the distance from the shower core, and t arr i va i is the arrival time of the considered 
secondary particle, as calculated by the code, and H is the injection point height in the 
atmosphere. The injection of the primary particle occurs by definition at zero time. 

We have recorded the delay distribution in different radial bins from the impact point 
of the shower axis for each different group of secondary particles. On the basis of the 
features of the lateral distribution of secondary particles, we use a logarithmic spaced 
radial zones, for a total number of 9 regions: 
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• I: 0-^4.64 m 



• II: 4.64^-10.00 m 

• III: 10.00^-21.54 m 

• IV: 21.54H-46.42 m 

• V: 46.42^-100.00 m 

• VI: 100. 00-^215.44 m 

• VII: 215.44-h464.16 m 

• VIII: 464.16^-1000.00 m 

• IX: 1000.00-h2154.42 m 

The contribution at distance exceeding the limit of region 9 are negligible for the 
considered energies. 

3.1 The arrival time probability distribution function 

We paid a particular attention to the choice of a suitable function which can be used to fit 
the delay distribution in any of the considered radial regions. If there is no interest in time 
delays larger than tens of nanoseconds, usually a good choice in experimental analysis is 
the gamma function, as in ref.|jj or |§j: 

r(i) = A ■ t p e- l - a (2) 

However, the authors often note how the gamma function cannot account for the long 
tail of the experimental delay distribution J2|. 

Our simulations also exhibit such long tails, but we have found that the use of a log- 
normal function allows a fairly satisfactory reproduction of the simulated data (at least 
when no smearing effects due to finite resolution of detector are included) in any of the 
aforementioned radial regions, up to a very long delay time (~ 10/is): 

/(t) = ^. e [~dog^) 2 /C] (3) 

The B and C parameters of such a distribution (essentially related to the mean delay 
and to its r.m.s) are also found to follow simple evolution as a function of shower energy 
and radial distance. 

The preference for the log-normal behaviour could reflect some intrinsic feature of 
the underlying processes. According to ref.|34j], the log-normal distribution arises when- 
ever we are in presence of a variable whose value takes a random proportion of that of 
the previous step in the stochastic process. However, we must stress that our choice of 
approximating function was an heuristic fact, aiming to the best approximation of the 
actual distribution. In fact we did not adopt any mathematical model as input for the 
shower development, thus we merely state that the log-normal distribution provides a 
better numerical approximation than the gamma function to the phenomenological delay 
distribution. Therefore, also the concept of goodness of fit has to be somewhat relaxed in 
the present discussion with respect to a rigorous statistical context. 
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Figure 1: Log-normal and Gamma fit to the delay distribution of secondary 7's and from 
primary 100 TeV protons, as detected at 2000 m a.s.l. in the V radial region (stand-alone 
FLUKA96 simulation. 



As an example, in Fig. |T| we show the time delay distribution for secondary gamma 
recorded at 2000 m a.s.l. in the radial region V, as produced by primary protons of energy 
100 TeV, fitted up to 200ns of delay, with log-normal and gamma functions. 

The differences are self evident and can be noticed at different energies, radial regions, 

etc. 

After this conclusion, we find more convenient to re-express all the results plotting log t 
instead of t and then fitting to a gaussian, as shown in the upper part of Fig. ^, which 
refers to the same case of Fig. ||. The distribution for secondary e + e~ has practically 
log-normal shape. It is also interesting to look at the distribution, for the same energy, 
radial bin, etc., when secondary muons and charged hadrons are selected. This is shown 
in the bottom part of Fig. [| 

The resulting distribution, for this primary energy and distance from the core, is even 
more asymmetric than that of electrons or photons. The consequence is that neither the 
pure gamma or the pure log-normal distribution succeed in fitting the delay distribution 
up to large t values. At smaller distances from the core, also for muons and charged 
hadrons, the log-normal approximation works better. However, in the following, we shall 
limit ourselves to the discussion of the results for e + e~ and 7's, since their density largely 
dominates that of muons, or residual charged hadrons. We recognize that the case of 
muons and other charged particle would deserve a dedicated study, also in view of the fact 
that it has been recently advocated the use of arrival time of muons for measurements 
related to the mass composition of primary cosmic rays, using the so called "Time Track 
Complementarity" (TTC) p5] . At present we are not able to perform with our tools a 
detailed study of this, since our Monte Carlo code does not allow, in this version, the 
treatment of nuclear projectiles and the primary energies required for this purposes. 
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Figure 2: Above: Gaussian fit to the log 10 t distribution of secondary 7's and from 
primary 100 TeV protons, as detected at 2000 m a.s.l. in the V radial region (stand- 
alone FLUKA96 simulation). Below: log 10 t distribution of secondary muons and charged 
hadrons, for the same primary energy and distance from the shower core. 



In any case, it already appears that also the log-normal fit is not sufficient to weight 
carefully the whole t range. Such a fit usually takes good care of the bulk of the events, 
but the extreme tails of the logi distributions exhibit deviations from a perfect gaussian 
behaviour. As a matter of fact, the left tail is very important, since it corresponds to the 
early front of the shower, that affects the experimental trigger time. The right tail refers 
to the very delayed component of EAS. In order to look for a more precise description 
of the simulation results, we characterise the log 10 1 distribution and its deviations from 
a pure log-normal one by means of a set of parameters. We find convenient the use of 
central moments (i.e. the moments calculated around the mean). Starting from the values 
of a finite set of these moments, the distribution can be reproduced ( and directly sampled 
in Monte Carlo) with a sufficient degree of approximation, as demonstrated in ref. J36[ |. 
We consider a variable £ defined by £ = igiioWzHi where ^\ and a are the average and 
the R.M.S. of the log 10 (i) distribution. If the /(£) distribution function (p.d.f. in the 
following) can be approximated by a standard normal <?(£) then /(£) can be expanded in 
series of the derivatives of the standard normal function: 

fiO = 0(0(1 + + (M4 - 3) H 4 + J^H 6 + ...} (4) 

Where H n are the Hermite polynomials of order n and \i n are the central moments of 

my- 

_ (logio (*i) - UlT 

i=l 

This is also known as Graham-Char lier expansion (GC in the following). We have 
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Figure 3: Log-normal and GC fits to the small delay region of the distribution of secondary 
7's and from primary 100 TeV protons, as detected at 2000 m a.s.l. in the V radial region 
(stand-alone FLUKA96 simulation. 



chosen to give the first four moments of log 10 (i) distributions. This set seems sufficient to 
express the features of the delay distributions in most cases, but there are some limitations, 
as discussed in the literature, which would demand for higher order moments. For instance 
when ^4 increases above a certain limit, moments of an order greater than n=4 must be 
considered in the expansion, otherwise the distribution becomes neither unimodal neither 
positive defined in the whole range. 

However we must stress how the moments of high order are subject to statistical 
fluctuations in case of a small number of entries in a distribution. In that case the errors 
induced by such fluctuations might be even larger than that due to the omission of such 
moments. 

To show how the G.C. expansion works, in Fig. ||] we show the comparison between 
the pure log-normal and the GC expansion fits to the short delay region of the same 
distribution of Fig. [l] (please note the linear scale!). Such an improvement can be seen 
also at large delays. 

In summary a better \ 2 1S obtained but usually the first two moments resulting from 
the fit, common to the two functions, are the same within the fit errors for the two 
case. We stress that, for most of the cases, these first two moments (i.e. the log-normal 
approximation) are sufficient to describe the time delay distribution from the experimental 
point of view (i.e. log(t/ns) >0), with the noticeable exception of long tails at large radial 
bins for proton showers, as will be discussed in Section 5. 

An important outcome of the series expansion in terms of moments is that it is possible 
to derive an expression which allows to have a direct Monte Carlo sample of the the desired 
p.d.f. with good accuracy starting from a random number normally distributed £ without 
any rejection. This is described and demonstrated in [p^| . The relation between variable 
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x to be generated with average /i, standard deviation a and central moments /%, and the 
normal distributed number £ is given by: 



x = \i + a 



£ + 



A*3 
6o- 3 



1 + 



24(7 4 



(« 3 - 3 f)-si«( 2 « 3 - 5 f 



(6) 



with the [i n giving the n-th central moments of the x distribution function. Of course \xi 



a 



An example of the functionality of this sampling method is shown in Fig. || where the 
M.C. output, the momentum expansion and the sampled data from this expansion are 
superimposed. 
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Figure 4: Comparison between the M.C. data the momentum expansion and the sampling 
described in 



Further differences, such as those discussed in Section || can be considered as an esti- 
mate of the systematic error related to the choice of different MC algorithms and interac- 
tion models. 



4 Results 

In this section we present the M.C. results giving information about the four moments that 
appear in the expansion of eq. |(| for different set of primary particle, secondary produced 
and distance from the shower core (radial bin). For practical convenience, instead of 
^2, /Mi in the following we shall use r.m.s (-s/JJa) , skewness (^3/^2 3//2 ) and kurtosis 
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(jL*4 / — 3) in place of n% and n±. For each set of the three quoted parameters, we obtained 
these quantities in two ways: by direct computation of the moments from the logarithm 
of the arrival times, and by fitting the log(t) distribution with the GC expansion with the 
moments used as fit parameters. A quantitative comparison of this two set of moments 
is made using the computed one as fixed parameters of the GC function and fitting the 
log(t) histogram with only the normalisation as a free parameter. In most cases, the two 
sets of resulting moments are in good agreement, within the errors, and also the \ 2 values 
from the two fits are similar. Important exceptions are those relative to the regions where 
we had low statistics. There, the moments obtained from the general fit generally give a 
worse approximation of the data than those directly computed. 

We will show how the behaviour of the parameters is rather smooth, so that all in- 
termediate cases in radius and detection height can be easily obtained by interpolation. 
We shall express the results in two ways: we will consider the moments of the average 
time distribution of single secondary particles (i.e. the distributions obtained summing 
the arrival times for each particles and for all the primary showers), but we shall also 
give quantitative information on the statistical fluctuations of these moments from event 
to event. Unfortunately, our simulation runs at the highest energies have only a small 
number of events, so that while the single particle distribution is always measured with 
high accuracy, the fluctuations on an event by event basis are studied with less precision. 

4.1 The p.d.f. moments 

The obtained moments]] for the distributions recorded at the observation level of 2000 m. 
a.s.l. are given in the tables reported in the Appendix. As previously stated, we limit 
ourselves to the results for secondary e + e~ and 7, which dominate the shower size. We 
refer mainly to the level of 2000 m a.s.l. since, in applications like gamma-astronomy, it is 
the most favourable among the three a-priori foreseen observation levels, although higher 
altitudes would be even more appealing for this purpose. The tables of parameters at the 
other lower observation levels can be obtained from the authors. 

In this section we present a subsample of these results in a graphical way in order 
to discuss the essential features. We plot the moment values without the corresponding 
error as derived from the fit. This is because those errors are dominated by the generated 
statistics and are lower than the actual shower to shower fluctuations, which instead are 
specifically reported in this paper. There is also another important comment about these 
results. Since the moments have been extracted from a fit, there exists some degree of 
correlation among them. This makes even less significant the error on a single parameter. 
We cannot report here the covariance matrix for all the relevant cases, and the irregularities 
visible in these figures are not simply due to statistical fluctuations. Correlations are 
such that the evolution of the resulting distribution is much smoother than that of single 
parameters. 

In Fig. [| the moment values are shown versus the logarithm of the primary energy for 
three radial bins, for the case of showers induced by 7 primaries and in Fig. |6| are plotted 
the same quantities for proton induced showers. 

The most striking feature is that the average value of log(i) is almost independent 
from primary energy (at least in the range considered by our simulation), kind of primary 
or secondary, and changes only with the radial bin. 

1 for the sake of simplicity, we shall use in the following the term moment, although r.m.s, skewness and 
kurtosis are not the true moments 
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Primary y -> y at 2000 m a.s.l. 
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Figure 5: First central momenta of the arrival time p.d.f. versus \og(E) for secondary 7 
from 7 primary in 3 different radial bins. Here r3, r5, and r8 stand for radial bin III, V 
and VIII. 
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If we adopt the log-gaussian approximation, we remind that, for this distribution, the 
correspondence between < t >, at and < log 10 (i) >, ^2 is given by: 

log 10 < t >=< log 10 (t) > +^ log e (10) (7) 

a t =< t > yj (10^ iog e (io) _ 1) ( 8 ) 

The increase in the average delay (which we remind is measured with respect to a light ray) 
as a function of distance from the shower core reflects the fact that moving to the outer 
regions of the shower, the average energy of secondary particles decreases, and therefore 
the deviations from a straight line trajectory become more important. 

The values of y/Ji2 are slightly different according to the different primary or secondary 
species, but is substantially constant with respect to the primary energy for 7 induced 
shower and shows a weak dependence from this parameter for proton primary. We also 
notice an apparent reduction of the variance (i.e. a shrinkage of the distribution) as 
function of the radial distance from the core. One should not be induced in error, since 
we plot the variance of logt: it can be verified that transforming back to the t variable, 
the effective width of the delay distribution (directly related to the time thickness of the 
shower disk) increases as a function of distance. Let us take a numerical example, for a 
primary gamma of 10 TeV, we have a r.m.s. around 0.6 in the third radial bin, when the 
average is 0.004. It means that the 68% fraction in logi is in the range 0.18 -j- 6.0 ns. 
Instead, in the eighth radial bin, the r.m.s. drops 0.28, while the average increases to 2.6. 
Therefore, the 68% fraction in logt is now in the range 118 -j- 1350 ns. 

The higher moments are very different for the two primaries: roughly speaking the 
proton induced showers have higher third and fourth moments, that is the log(i) p.d.f. 
has a less log-normal shape. In order to understand that, but also to give a clearer idea 
of the evolution of time delay distributions, we show in Fig. the log t distributions at 
10 TeV for 7 and proton showers in few different radial bins. We can clearly see that 
while there is a nearly log-gaussian bulk almost identical in the two cases, an additional 
tail at high delays appears for proton showers, particularly visible at large distance from 
the shower core. Our fit is unable to follow completely this tail, which presents itself as 
a distinct family of particles. Indeed the population of these highly delayed particles is 
only a very small fraction of the total number of arriving e.m. particles. We postpone the 
specific discussion on this phenomenology in section ||. 

Coming back to the general features of the delay distribution, we stress that the ex- 
tracted moments describe the average arrival time distribution functions. Obviously, re- 
markable fluctuations can be detected in a event-by-event analysis of the generated data. 
For example, in Fig. ^ we superimpose the electron arrival time distribution for 4 different 
showers induced by a 100 Tev 7 primary. The differences in shape and normalisation 
for the different showers are evident. This effect is amplified in the showers induced by 
proton primaries. This is highly correlated to the large fluctuations in the height of the 
first inelastic interaction. In Fig. ^ are presented the momenta fluctuation for 7 primaries, 
to be compared with much bigger fluctuation presented in Fig. [K^, corresponding to the 
proton primaries. As expected, all the fluctuations shows a decreasing behaviour with the 
energy. 
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Primary Protons -> yat 2000 m a.s.l. 
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Figure 6: First central momenta of the arrival time p.d.f. versus log(-E) for secondary 7 
from proton primary in 3 different radial bins. Here r3, r5, and r8 stand for radial bin III, 
V and VIII. 
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Figure 7: Logiot distributions at 10 TeV for proton (above) and 7 (below) showers in few 
different radial bins. Here r2, r5, and r8 stand for radial bin III, V and VIII. 
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Figure 8: Electrons arrival time distributions for single showers from gamma 100 TeV 
primary. 
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Primary y -> y at 2000 m a.s.l. 
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Figure 9: Fluctuations of first central momenta of the arrival time p.d.f. versus log(-E) 
for secondary 7 from gamma primary. Here r3, r5, and r8 stand for radial bin III, V and 
VIII. 
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Primary Protons -> yat 2000 m a.s.l. 
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Figure 10: Fluctuations of first central momenta of the arrival time p.d.f. versus \og(E) 
for secondary 7 from proton primaries. Here r3, r5, and r8 stand for radial bin III, V and 
VIII. 
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Rad. bin time thr. N ^ l ^ d x 10" 4 M 7 N„ eutrona n K+e - 

-tv total 1* delay ed ^delayed ^delayed 

H > 10' 2 ns 9^3 89.5 (%) 3M%) 7.0(%) 
V > 10 3 ns 5.4 30.0 (%) 59.0 (%) 4.0(%) 
VIII > 10 4 ns 196.8 6.5 (%) 93.0 (%) 0.5(%) 

Table 1: Percentage of delayed particle and relative abundance of neutrons, photons and 
electrons for different radial bins, in EAS induced by a 100 TeV proton primary. 



5 Discussion of the origin of delayed particles 

The time structure of e + e~ and 7 secondaries from hadronic and pure e.m. showers shown 
in Fig. ^ is rather similar, as far as the nearly log-gaussian bulk is concerned. The main 
difference consists in the appearance, in the hadronic showers, of a slow e.m. component 
which becomes more and more visible at large distance from the shower core. Such a 
component seems to constitute a different family, really well separated from the major 
population. It can extend itself up to tens of microseconds. We have investigated in some 
detail the simulation steps, in order to understand the physical origin of such a delayed 
component. We have learned, of course, that the bulk of our time delay distribution 
come from the relativistic part of the shower, namely from the e.m. showers coming from 
ttq decay. Instead, the delayed family in hadronic showers comes from the de-excitation 
of air nuclei after the interactions with neutron of energy below few tens of MeV. These 
neutrons are quite slow, and for this reason separate easily from the relativistic component. 
Of course, this secondary electro-magnetic component is a small fraction of a much larger 
population of neutrons which dominates the shower disc at large distances. Table |l| shows 
the percentage fraction of delayed particles, neutrons, 7 and electrons, above a radial 
dependent time threshold, for 100 TeV proton primary at 2000 m a.s.l.. The time threshold 
for a given radial bin has been chosen in such a way to approximately define the transition 
region between the log-normal bulk and the additional tail. 

The treatment of nuclear effects, and in particular of de-excitation processes, is one 
of the most interesting feature offered by the FLUKA code. The description of the mod- 



elization of these processes is given in ref. [37] and |p8| . We have to notice that also in 



the FLUKA interface of GEANT 3.21 a nuclear evaporation model is present, but it much 
more simplified and approximated with respect to that of the quoted references. We are 
aware that this new evaporation model is now also interfaced to CORSIKAf^], version 
5.60, when the DPMJET interaction model|39| is chosen. However such a package in the 
frame of CORSIKA is not yet available for public distribution. 

Here we try to summarize the possible different low-energy reactions which can con- 
tribute to the process under investigation. 

a) capture reactions (n, 7): very probably these are not the relevant phenomena in our 
case, since these would give photons with much larger delays than those considered 
in our distributions. However, in general they have a small cross section. This 
reaction can be relevant for Ar (630 mbarn) through the process 40 Ar(n, 7) 41 Ar. On 
this respect we have to make the following remark. The present simulation has been 
performed in the case of a dry atmosphere, but in case of a relevant presence of 
humidity in air or near the soil (in presence of snow for instance), the build up of 
capture reactions can be relevant due to the thermalization of neutrons on H nuclei. 
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b) reactions (n, n'), where the neutron is scattered leaving the nucleus in an excited 
state. The residual nucleus will decay to the ground state emitting a cascade of 
photons between the different levels. This phenomenon is one of the most important 
sources of delayed photons, since an energy just above the first excited level is enough. 
For N, the first three excited levels are at 2.319, 3.948, 4.915 MeV. For Ar they are at 
1.461, 2.121, 2.524 MeV, while for O they are at 6.049, 6.130, 6.919 MeV (therefore 
is less relevant in the O case). The excitation of higher levels is usually negligible. 
The lifetime of these excited levels is quite short (ns scale or less), so that these 
photons are promptly emitted. Furthermore, the model makes use of an isotropic 
angular distribution. This is not very different from reality. 

c) There are also reactions with production of other particles together with additional 
7's. We can quote, in order of increasing energy threshold: 

(n,p); among these, an important exception is the celebrated 14 N(n,p) 14 C, which 
does not produce 7's; 

(n, 2n). This is practically closed for O nuclei below about 20 MeV, but is open for 
the other nuclei of air even at lower energies; 

(n,a); 

(n,d); 

(n,t); 

(n,kp), with k produced protons, 
(n,ka), with k produced a's 

For each of these processes, the cross section rises quite fast above the threshold. 

The reason why we notice the delayed family at large distance is that the low energy 
neutrons can be found more easily in the periphery of the shower, where energy is degraded. 
Furthermore, this non relativistic component is obviously much more smeared in space 
than the relativistic one. Of course, the delayed neutron component can give much larger 
signals in a given detector (for instance in a scintillator array) with respect to the associated 
secondary e.m. component. However, we would like to stress that in case of nuclear 
projectiles (a case outside the scope of the present simulations) we shall have also excitation 
of the nuclear fragments coming from the projectile. Unlike the case of the excitation of 
the target nuclei, the photons (and other particles) following the de-excitations of the fast 
moving fragments will be Lorentz-boosted in the laboratory frame, so that to contribute 
in a more efficient way to the energy deposition. In particular we advance the hypothesis 
that they can give rise to subshowers having some delay with respect to the first particles 
in the EAS disc. For instance, we suggest that this could be the basis for an explanation 
of the observation of a delayed component |24|] in VHE Extensive Air Showers, as observed 
in the COVER-PLASTEX detector in the GREX array |]. 

We wish to point out that the inclusion of these processes in simulation codes for 
high energy cosmic ray physics has not yet become a common practice. Therefore, these 
phenomena might have escaped from other simulation studies, while the existence of such 
delayed particles might have some experimental relevance. 
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6 Comparison with experiments 



We compare our results with the EAS-TOP |4j array experimental data. In Fig. 11 we 
plot the average secondary delay from the first particle, with respect to the radial bin. 
This values refer to the photons produced by a 100 TeV, proton initiated EAS as detected 
at 2000 m a.s.l. The EAS-TOP signal is dominated by the secondary electrons, but we 
remark that, at least in this range of parameters, the differences with respect to the time 
distribution of secondary electrons appear to be very small. On the other hand the number 
of photons is much larger^], and this allows a reduced statistical uncertainty. The average 
value we obtain seems to reproduce the EAS-TOP experimental data. These are divided 
in 2 samples, corresponding to the cases in which they had at least 1 or 4 particles in 
each of the triggered sub-detectors. The event class corresponding to at least 1 particle 
in the detector should correspond to a primary energy close to the trigger threshold of 
EAS-TOP (~ 100 TeV, mostly from proton primary). At distances from the core lower 
than 100 m, the two classes are not distinguishable. In the same figure we also show that 
our results are systematically higher than those of ref. [|J . There the calculation was done 
only for gamma primaries, but our results, in the same range of distance from the shower 
core, do not exhibit substantial differences between the two species of primaries. Also 
the delay distribution width (that is the shower disk thickness), is well reproduced by our 
simulation: in Fig. |ll| is reported the width of the delay distribution as measured at 70% 
of the height of the peak, for different radial bins. The agreement between Monte Carlo 
and the EAS-TOP data is very good. In spite of the goodness of the simulation results 
we want to stress the effect of a finite detector resolution on the simulated arrival time of 
the secondaries. In particular, in the first radial bins (near the axis of the shower), the 
arrival time distribution are peaked at time values which can be lower than the typical 
experimental time resolution. Thus the result of the folding of the simulated p.d.f. with 
the experimental error distribution drastically modifies the distribution in these radial 
bins, producing a flattening at arrival times close to the time resolution value. 



This effect is shown in Fig. |12J where we plotted a log-normal function with < 



log 10 (t/ns) >= —0.54 and 0"z O9lo ( t / ns ) = 0.75 (similar to the delay distribution expected 
at 2000 m a.s.l for secondary photons from primary gamma, in the 2nd radial bin), super- 
imposed to the same distribution folded with a gaussian error function with a exp = 1 ns. 
The effect of the experimental resolution is clearly seen on the bulk of the events, even if 
the average value of the arrival time is bigger then the resolution value. 

We remark also that usually the experimental results are expressed as particle delay 
distribution as measured since the arrival of the experimental trigger. This is different 
from the delay time from the light cone that we used in M.C. These two distributions are 
almost equal if the shower is so populated that the first triggering particle has a minimum 
delay from the light cone, but can be quite different for low statistic showers, where the 
triggering particle may arrive after a sizeable delay from the light cone. In this case the 
M.C. distribution are biased toward the longer delays with respect to the experimental 
results by an amount of time given by the average delay of the first particle. 



7 Conclusion 

We have simulated the time delay of vertical Extensive Air Showers by means of the 
FLUKA Monte Carlo code. We have shown how the logarithm of the time delay dis- 
tribution of a shower front can be easily expressed as a gaussian expansion in terms of 
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a few parameters. Such parameters exhibit a smooth behaviour as a function of energy, 
thus allowing the extrapolation to higher energies of the presented Monte Carlo simulation 
results. 

The proposed expansion allows the construction of fast algorithms useful in detector 
design and data interpretation, when the running of a complete detailed simulation would 
result in unnecessary heavy job. Differences between protons and gamma initiated show- 
ers are observed. In particular, hadronic showers exhibit a delayed component at large 
distance from the shower axis. These differences have been interpreted as due to effects 
coming from non relativistic particles in the shower. In particular, we found that the de- 
excitation of air nuclei following the interaction with low energy neutrons are the essential 
contributions to the highly delayed e.m. component in hadronic showers. We remark how 
this can be obtained at simulation level, only with very specialized and detailed codes, not 
commonly used so far in cosmic ray physics. In particular situations, these results can be 
used as an interpretation basis for the understanding of apparently anomalous delays. 
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Appendix: Tables of Moments and of their Fluctuations 

Table 2: Central momenta for the distribution function of log (^) for secondary 7 pro- 
duced at 2000 m a.s.l. by a 7 primary, as obtained from the FLUKA stand-alone simulation 



Primary 7, Sec. 7, 1^=2000 m a.s.l. 



E = 562. GeV 



Rad. 


bin 


average 


sigma 


skewness 


kurtosis 


1 




-.77± .45 


1.03± .28 


.3± .7 


3.2± 1.1 


2 




-.35± .30 


.83± .19 


.5± .6 


3.0± .7 


3 




.05± .26 


.71± .17 


.7± .7 


4.1± 1.3 
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.51± .25 


.61± .15 
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3.1± .9 
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E 
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bin 
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skewness 


kurtosis 
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3.1± 1.0 
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sigma 
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Table 3: Central momenta for the distribution function of log (^) for secondary e + e~ 
produced at 2000 m a.s.l. by a 7 primary, as obtained from the FLUKA stand-alone 
simulation 



Primary 7, Sec. e + e , h^ et =2000 m a.s.l. 
E n = 562. GeV 
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Table 4: Central momenta for the distribution function of log (^) for secondary 7 pro- 
duced at 2000 m a.s.l. by a proton primary, as obtained from the FLUKA stand-alone 
simulation 



Primary protons, Sec. 7, 11^=2000 m a.s.l. 
E n = 562. GeV 
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.2± .1 


9 84- 3 

Z.Q4Z .O 


5 


1.06± .03 


.49± .02 


.13± .06 


2.7± .3 


6 


1.57± .07 


.43± .02 


-.2± .2 


2.7± .5 


7 


2.06± .05 


.37± .02 


-.2± .3 


2.8± 1.0 


8 


2.51± .06 


.33± .03 


-.3± .4 


3.4± 1.4 


9 


2.97± .03 


.44± .03 


.3± .5 


5.9± 1.2 




Eo 


= 31620. GeV 




Rad. bin 


average 


sigma 


skewness 


kurtosis 


1 


-1.34± .52 


.92± .07 


.4± .4 


3.3± .5 


2 


-.56± .24 


.74± .08 


.5± .2 


3.6± .5 
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3 -.02± .12 .64± .06 .3± .2 3.3± .3 

4 .53± .04 .56± .04 .25± .07 3.0± .2 

5 1.07± .05 .49± .03 .0± .1 2.7± .3 

6 1.59± .05 -.43± .02 .1± .2 2.7± .4 

7 2.09± .05 .37± .03 -.3± .3 3.0± 1.1 

8 2.53± .06 .33± .04 -.3± .5 3.4± 1.2 

9 3.00± .08 .42± .03 .7± .4 5.5± .8 

E = 56230. GeV 

Rad. bin average sigma skewness kurtosis 

1 -1.29± .33 .88± .05 .5± .2 3.4± .3 

2 -.58± .05 .70± .05 .51± .09 3.5± .2 

3 -.03± .03 .62± .04 .46± .03 3.3± .2 

4 .52± .04 .55± .02 .27± .06 3.0± .1 

5 1.08± .03 .48± .02 .12± .05 2.7± .1 

6 1.61± .03 .41± .01 -.1± .1 2.7± .3 

7 2.10± .04 .36± .02 -.3± .1 2.9± .6 

8 2.56± .04 .31± .02 -.2± .4 3.2± .8 

9 3.05± .05 .43± .02 .7± .2 5.8± .5 

E = 100000. GeV 

Rad. bin average sigma skewness kurtosis 

1 -1.25± .11 .88± .05 .5± .1 3.3± .2 

2 -.58± .08 .71± .06 .5± .1 3.7± .2 

3 -.03± .04 .64± .03 .39± .09 3.4± .2 

4 .52± .02 .55± .02 .28± .05 3.04± .07 

5 1.08± .03 .48± .02 .08± .02 2.8± .1 

6 1.61± .03 -.41± .01 .11± .05 2.7± .2 

7 2.10± .02 .35± .01 -.3± .1 2.9± .6 

8 2.56± .03 .31± .02 -.2± .2 3.2± .6 

9 3.02± .04 .40± .02 .6± .2 5.1± .6 
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Table 5: Central momenta for the distribution function of log (-M for secondary e + e~ 
produced at 2000 m a.s.l. by a proton primary, as obtained from the FLUKA stand-alone 
simulation 



Primary protons, Sec. e + e , hrf et =2000 m a.s.l. 
E n = 562. GeV 



Rad. bin 


average 


sigma 


skewness 


kurtosis 


1 


.26± .30 


.43± 


.10 


.37± .7 


1.90± .10 


2 


.49± .27 


.63± 


.11 


.85± .6 


3.74± .8 


3 


.57± .27 


.59± 


.11 


.2± .6 


2.2± .8 


4 


.66± .26 


.54± 


.11 


.2± .7 


2.1± 1.4 


5 


.93± .28 


.50± 


.12 


-.2± .8 


3.4± 2.4 


6 


1.30± .23 


.43± 


.12 


.0± .7 


2.9± 1.7 


7 


1.72± .28 


.43± 


.19 


.1± 1.0 


3.5± 1.8 


8 


2.23± .36 


.44± 


.20 


.4± 1.0 


3.5± 1.5 


9 


2.74± .34 


.60± 


.27 


.3± 1.0 


5.8± 5.8 




E 


= 1000. GeV 




Rad. bin 


average 


sigma 


skewness 


kurtosis 


1 


.17± .29 


.69± 


.19 


.2± .5 


2.2± .6 


2 


-.03± .36 


.79± 


.15 


.2± .5 


2.1± .8 


3 


.17± .25 


.68± 


.11 


.2± .5 


2.3± .6 


4 


.60± .17 


.54± 


.11 


-.2± .6 


2.4± .8 


5 


1.02± .26 


.47± 


.12 


-.1± .6 


3.3± 1.0 


6 


1.34± .24 


.43± 


.12 


-.1± .8 


3.6± 3.1 


7 


1.81± .28 


.58± 


.17 


-.2± 1.1 


5.8± 3.7 


8 


2.32± .34 


.56± 


.19 


1.0± 1.0 


5.9± 2.1 


9 


2.77± .37 


.50± 


.18 


.2± .9 


3.2± 1.9 




Eo 


= 1778. GeV 




Rad. bin 


average 


sigma 


skewness 


kurtosis 


1 


.42± .40 


.62± 


.20 


.54± .5 


2.50± .4 


2 


.15± .40 


.78± 


.19 


-.2± .4 


2.3± .3 


3 


.28± .21 


.69± 


.13 


.0± .5 


2.4± .7 


4 


.56± .18 


.56± 


.13 


-.1± .6 


2.0± .8 


5 


.95± .18 


.46± 


.11 


.0± .7 


2.8± 1.6 


6 


1.33± .19 


.41± 


.10 


.2± .7 


3.2± 1.9 


7 


1.81± .20 


.36± 


.11 


.3± .8 


3.1± 1.9 


8 


2.31± .29 


.51± 


.16 


.2± 1.0 


4.8± 2.3 


9 


2.86± .24 


.55± 


.20 


1.1± .9 


5.1± 1.8 




Eo 


= 3162. GeV 




Rad. bin 


average 


sigma 


skewness 


kurtosis 


1 


-.32± .30 


.77± 


.20 


.1± .5 


2.2± .8 


2 


-.14± .31 


.80± 


.14 


-.2± .4 


2.7± .9 


3 


.20± .23 


.67± 


.11 


-.2± .4 


2.5± .5 


4 


.63± .18 


.53± .09 


-.1± .4 


2.8± 1.0 
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E 





i m 4- is 

l.Ul± .15 


.47± .11 


-.2± .4 


9 Q4_ ft 

z.y± .d 


ft 





1 qft-L- 1 E 
l.O0± .10 


.40± .09 


-.2± .8 


q Q4- q Q 

o.o± o.o 


7 




1 894- 1 Q 

i.oz± .iy 


.37± .13 


-.1± 1.0 


q 1 4_ q « 
o.i± o.o 


o 
o 




O QO-L Oft 
Z.OZ± .ZD 


.42± .18 


.3± .8 


q 74- O 1 
O. ( ± z.l 


n 




884- Q7 

z.oo± .0/ 


.58± .20 


.0± .7 


E E I 1 E 
0.O± 1.0 






E 


= 5623. GeV 




Kad. 


bin 


average 


sigma 


skewness 


kurtosis 


i 
i 




91 _|_ 97 
-.Ollt .Zl 


.85± .20 


-.1± .6 


9 /i4_ i n 

Z.4zt l.U 


9 
z 




9D4- zLfi 


.86± .16 


-.1± .4 


9 ^4- 
z.ozn .u 


q 
o 




qi _l 9/1 
.ol± .Z4 


.63± .14 


-.1± .6 


9 /I4- 1 1 
Z.4± 1.1 


1 

4 




E.Q4- 9ft 

.oy± .ZD 


.55± .11 


-.4± .4 


9 04- E 

z.y± .0 


E 






i nn4- i /i 

l.UU± .14 


.47± .10 


-.2± .5 


9 84- 9 1 
Z.o± Z.l 


ft 





1 OO-L 1 E 

l.zy± .lo 


.50± .11 


-.5± .8 


/I 74_ O A 
4. ( ± Z.4 


7 




1 8q4- 1 Q 

l.oo± .lo 


.40± .11 


.0± .5 


q 84- i e 
o.o± 1.0 




o 




9 q/14- 9n 

Z.o4± .zu 


.41± .16 


.4± .9 


A 1 4- 9 1 
4.1± Z.l 


n 




Z.OO± .DO 


.65± .19 


.5± .9 


/l q_i_ i ft 
4.o± 1.0 






Eo 


= 10000. GeV 




Kad. 


bin 


average 


sigma 


skewness 


kurtosis 


i 
i 




i K4- q9 

-.10± .oZ 


.83± .19 


-.3± .4 


9 74- 7 
Z. f ± . 1 


o 
z 




1 E_|_ 9Q 

-.loin .zo 


.79± .15 


-.2± .4 


1 84- E 
l.ozn .0 


q 

o 




ooi no 

.zo± .zz 


.68± .10 


-.3± .3 


o q_i_ ft 
z.o± .0 


/I 




ftq_l_ 1 O 
.0o± .lZ 


.55± .08 


-.2± .3 


O 74- q 

z. ( ± .6 


E 






l.Ul± .1U 


.46± .06 


-.2± .5 


9 Q4_ 1 /I 

z.y± 1.4 


ft 





1 q84- 1 Q 

l.oo± .lo 


.39± .07 


-.1± .8 


q 94- A A 
o.Z± 4.4 


1 




i qi i in 
l.olrt .1U 


.37± .07 


.1± .9 


q q_L /i i 
o.o± 4.1 


g 

o 




O QQ4- 1 A 
Z.OO± .14 


.38± .09 


.2± .9 


o 1 Oft 

0.4± Z.O 


Q 

y 




9 884- 99 
z.oo± .ZZ 


.56± .15 


.5± .7 


E. 74- 1 ft 
0. 1 ± 1.0 






Eo 


= 17780. GeV 




i>] 

Kad. 


bin 


average 


sigma 


skewness 


kurtosis 


i 
i 




994_ /|9 
-.ZZzt .4Z 


.82± .14 


-.1± .3 


9 74- 
Z. 1 zt .o 


9 




1 94. 9« 
-. izzn .zu 


.78± .08 


-.2± .5 


44- 1 1 
o.^tzn i . l 


q 
o 




oi4_ no 
.zl± .uy 


.67± .06 


-.1± .2 


O 04- A 

z.y± .4 


4 






.53± .05 


-.24± .10 


2 7+3 


5 




1.06± .04 


.43± .04 


-.1± .2 


2.8± .5 


6 




1.40± .06 


.37± .03 


.0± .5 


2.9± 3.5 


7 




1.85± .06 


.34± .03 


.1± 1.3 


3.4± 5.8 


8 




2.34± .14 


.38± .09 


.1± .8 


3.8± 2.1 


9 




2.93± .13 


.50± .08 


1.1± .5 


4.5± 1.5 






Eo 


= 31620. GeV 




Rad. 


bin 


average 


sigma 


skewness 


kurtosis 


1 




-.78± .43 


.91± .13 


.0± .3 


2.0± .5 


2 




-.14± .21 


.71± .08 


-.1± .2 


2.6± .3 
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o 
o 


004- 1 

.zy± .iz 


R94- OR 

.oz± .uo 


04- Q 
-.Z± .o 


O 84- 8 
Z.O± .8 


A 

4 


7n4- r\K 

. iu± .uo 


.0Z± .Uo 


i A-L 1 n 
-.14± .1U 


O 74- Q 
Z. 1 ± .O 


K 




1 D7-I- Dfi 
l.UlIt .UO 


AA-i- n^ 

.44zt .UO 


94- 1 
-.Zzt .1 


o n4- fi 

O.Uzt .0 





1.4o± .UO 


Q74- HQ 

.oi± .Uo 


14- Q 
-.1± .O 


Q n_L 1 o 
o.U± l.Z 


7 


1 874- n^ 
1.0/ ± .UO 


.o0± .Uo 


.u± .y 


Q 04- Q A 
O.Z± 0.4 


Q 

o 


Z.oo± .Uo 


QO l H7 
.oo± .III 


Q I Q 

.o± .o 


/( 04- O Q 
4.Z± Z.o 


n 

y 


Z.OO± .10 


/1 7_1_ HQ 

.4/± .Uo 


FT I r 
.0± .0 


0.0± l.U 




E 


= 56230. i 


GeV 




Rad. bin 


average 


sigma 


skewness 


kurtosis 


i 
i 


7Q4- OQ 
-. 1 OzSz .ZO 


Qn_i_ nfi 
.yuzt .uo 


14- 9 
.lit .z 


1 04- 

i.yzt .o 


9 


-. ± ( nz .u^ 






9 ^4- ^ 
z.ozn .o 


Q 

o 


oox n^ 

.zy± .uo 


fil 4- HQ 

.01± .Uo 


no4- n7 
-.UZ± .U/ 


o oo4- no 
z.zy± .uy 


4 


.oy± .U4 


_i_ no 
.01± .UZ 


1 c_i_ no 
-.lo± .Uo 


O A 4_ /I 

Z.4± .4 


K 




i ns4- nQ 
l.Uo± .Uo 


a Q-i- no 

.4o± .UZ 


14- 1 
-.1± .1 


O 74- /I 
Z. 1 ± .4 





1 /1£4- HQ 

1.40± .Uo 


oc i no 
.o0± .UZ 


n i q 
.U± .o 


Q 04- 1 o 

o.u± i.y 


7 


1 094- O/l 

i.yz± .U4 


Q/i 4- no 

.o4± .UZ 


94- ^ 
.Z± .O 


Q Q4- 9 1 

o.o± z.i 


g 
o 


ZA l ii± .UO 


A A-L C\A 
.44± .U4 


Q 1 O 

.o± .Z 


/i 04- 1 n 
4.y± i.u 


o 

y 


o o24- no 

z.yo± .uy 


ci i n/i 

.01± .U4 


^4- Q 
.0± .O 


0.1± .0 




Eo 


= 100000. 


GeV 




Rad. bin 


average 


sigma 


skewness 


kurtosis 


1 


-.72± .14 


.90± .03 


.1± .2 


1.9± .2 


2 


-.14± .08 


.70± .06 


.0± .1 


2.3± .1 


3 


.29± .04 


.61± .03 


-.13± .07 


2.6± .1 


4 


.69± .02 


.52± .02 


-.14± .04 


2.68± .09 


5 


1.08± .03 


.43± .01 


-.08± .06 


2.8± .2 


6 


1.45± .03 


.36± .01 


.0± .2 


3.0± 1.9 


7 


1.90± .03 


.34± .02 


.0± .4 


3.2± 2.1 


8 


2.42± .04 


.37± .03 


.4± .2 


4.0± 1.2 


9 


2.96± .07 


.51± .05 


.6± .3 


5.7± 1.1 
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